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Using a thermodynamic integration scheme, we compute the free energy cost per unit area, 7, 
of forming an interface between a crystal and a frozen structured wall, formed by particles frozen 
into the same equilibrium structure as the crystal. Even though the structure and potential energy 
of the crystalline phase in the vicinity of the wall is same as in the bulk, 7 has a non-zero value 
and increases with increasing density of the crystal and the wall. Investigating the effect of several 
interaction potentials between the particles, we observe a positive 7 at all crystalline densities if the 
potential is purely repulsive. For models with attractive interactions, such as the Lennard-Jones 
potential, a negative value for 7 is obtained at low densities. 


I. INTRODUCTION 

The excess free energy of a crystal in contact with a 
wall is an important quantity governing many interfa¬ 
cial phenomena such as wetting and heterogeneous nu- 
cleation [H-Q. Theoretical approaches based on cell the¬ 
ory have been used to compute this quantity for model 
systems such as hard spheres Q. Atomistic simulation 
studies have also been carried out employing widely-used 
interaction potentials such as Lennard-Jones (LJ) models 
and systems of hard spheres iOl. 

In most previous studies, the cry^al is either in contact 
with a flat structureless wall or a structured 

wall consisting of particles attached to sites correspond¬ 
ing to an ideal crystalline structure Q- In both of these 
situations, the structure and energy of the crystalline lay¬ 
ers in the neighborhood of the wall is not the same as in 
the bulk leading to an excess free energy as compared to 
the bulk. 

An interesting situation arises when the energy and 
structure of the crystalline layers in contact with the wall 
remains the same as in the bulk. Such a scenario occurs 
when the wall consists of particles frozen into positions 
occupied by the same crystalline phase at equilibrium. 
In this case, the wall is a frozen version of the same crys¬ 
tal. Then, it is of interest from a fundamental statistical 
mechanical point of view, to determine the excess free 
energy of the system in comparison to the bulk. 

An analogous situation corresponding to a supercooled 
liquid in contact with a wall consisting of particles frozen 
into the same amorphous structure as the liquid has been 
widely studied in recent years in order to investigate 
the relaxation dynamics of glass-forming systems | 12 - ll^ . 
Recently, an experimental investigation of such a system 
was carried out using colloidal particles [ 2 ^. Such com¬ 
putational and experimental studies extract a growing 
length scale based on a slowing down of the dynamics 
near the wall and relate it to the slowing down of glassy 
dynamics in the bulk. In these investigations, it is implic¬ 
itly assumed that the thermodynamics of the system is 
unperturbed by the wall if an average is carried out over 
the thermal fluctuations and different realizations of the 
wall However, in a recent work we obtained a 


non-zero interfacial free energy for a glass-forming binary 
Lennard-Jones liquid (2lj | in contact with a wall with the 
same amorphous structure as the supercooled liquid, in¬ 
dicating that thermodynamics of the liquid is affected by 
the presence of such a frozen wall . 

In this work, we compute the excess free energy of a 
crystal in contact with a wall formed by its frozen version, 
using molecular dynamics simulation s 12311 in conjunc¬ 
tion with thermodynamic integration [24| . We obtained 
a non-zero interfacial free ener gy f or particles interacting 
via a Lennard-Jones potential |25l |. which increases with 
increasing density of the crystal suggesting that the ther¬ 
modynamics of the crystal in contact with such a wall is 
perturbed and not the same as in the bulk. Moreover, 
at low densities 7 is negative, indicating an increase of 
entropy as compared to the bulk. To test the robustness 
of our results and find out whether a negative value of 7 
is model-dependent or not, we carried out investigations 
with two other interaction potentials: a purely repulsive 
inverse-twelfth power soft-sphere potential and the hard 
sphere interaction (using a potential) [ 2 ^. How¬ 

ever, 7 for the two purely repulsive interactions, while 
non-zero, is positive at all crystal densities. 

In the next section we specify the interaction poten¬ 
tials considered in this work. In section imi the system 
is described and then we outline the thermodynamic in¬ 
tegration scheme adopted to compute 7 . In section El 
we describe the simulation details and then present our 
results in section IVTl Finally, we end with a conclusion. 


II. INTERACTION POTENTIALS 

In order to study the robustness of our findings, we ob¬ 
tained 7 corresponding to three different interaction po¬ 
tentials. The interaction potentials are denoted by U{(r) 
[i = 1,2,3], r being the distance between the particles. 
Energies and lengths are expressed in units of the param¬ 
eters e and a, respectively (see below), and the masses 
of the particles are set to m = I. In the following, all 
other quantities are given in terms of the latter param¬ 
eters, i.e. time t in units of r = pressure P in 

units of e/(T^, and the interfacial free energy in units of 
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The following interaction potentials are considered in 
this work: 

(i) A modified Lennard-Jones interaction potential 
(inLJ) is chosen, defined by 


ui(ry) = (t){nj) + Cl 


with 


(j){rij) =4e 




( 1 ) 

( 2 ) 


for 0 < Vij < 2 . 3 ( 7 , and 

uiir^j) = C 2 +C 3 +C 4 + C 5 (3) 

for 2 . 3(7 < rij < = 2.5a and u{rij) — 0 for > 

^cut_ constants in Eqs. © and m are given by 

Cl = 0.016132 e, C 2 = 3136.6 e, C 3 = -68.069 e, C 4 = 
-0.083312 e, and C 5 = 0.74689 e. 

(ii) The second potential we consider is the purely 
repulsive force-shifted inverse-twelfth power potential 
(fsil 2 ), 


U 2 {rij) = ipinj) - (/>(rc) - (r 


d^pinj) 


dnj 


( 4 ) 


where, V'(nj) = 4 . 0 e(c 7 /rij)^^. Here, the potential is trun¬ 
cated and shifted such that U 2 {rij) and its derivative with 
respect to goes to zero at rc = 2.5a. 

(iii) The third potential is a hard-sphere interaction 
modeled by a inverse-power potential (1256): 

U3inj) = e{ —) . (5) 

The potential is cut off at rc = 1.2a. As shown previously 
[ 2 ^, the coexistence and other bulk properties of this 
interaction potential are very similar to those of the hard- 
sphere system. 

In Fig. [U we show the three interaction potentials as 
a function of the distance r between the particles. 


III. SYSTEM SETUP 

Our systems consist of face centered cubic (fee) crys¬ 
tals of N particles, with the (100) face oriented along 
the z direction, enclosed in a simulation box of dimen¬ 
sion Lx X Ly X Lz. Periodic boundary conditions are 
maintained along the x and y directions. Along the 
z direction, the particles are confined by several layers 
of the same crystalline particles frozen into their equi¬ 
librium configuration (see Fig. HI). Since there are two 
wall-crystal interfaces, the total area of the interface is 
A = 2LxLy. In order to prevent any crystalline particle 



FIG. 1. (Color online) Interaction potentials corresponding 
to a modified LJ potential (mLJ), the inverse-256th power 
interaction approximating the hard sphere potential (1256), 
and the force-shifted inverse-twelfth power potential (fsil 2 ). 


from penetrating the wall, a short-ranged flat wall, mod¬ 
eled by a Gaussian potential is inserted at the boundaries 
along the z direction (at z = 0 and Lz). As shown in our 
earlier works [ 2 ^, [ 2 ^ , such a short-ranged flat wall leads 
to a negligible contribution to the interfacial free energy 
(r; O.OOI 7 ). 


IV. THERMODYNAMIC INTEGRATION 

The interfacial free energy is defined as 7 = (Lgystem — 
A’buik)/A, where Lljystem is the free energy of the crystal 
confined by the frozen walls on both sides, while Lbuik is 
the free-energy of the bulk crystal with periodic bound¬ 
ary conditions along all the three Cartesian axes. Our 
objective is to compute 7 using molecular dynamics sim¬ 
ulation and thermodynamic integration. 

The thermodynamic integration scheme adopted in 
this work is based on a similar approach devised by us 
in previous studies to compute the interfacial free energy 
of a liquid/crystal in contact with structured walls [ 3 -Q 
and to determine the crystal-liquid interfacial free en¬ 
ergy m, [ 2 ^. The goal is to transform a bulk crystal 
with periodic boundaries in all three Cartesian-axes di¬ 
rections into a crystal in contact with frozen crystalline 
layers comprising the frozen walls on either side. To ac¬ 
complish this, it is neccessary that the boundary condi¬ 
tions along the z direction are rearranged at some point 
during the transformation and that the crystalline parti¬ 
cles do not cross the boundaries during this rearrange¬ 
ment. Hence, in the first step, two extremely short- 
ranged flat walls, modeled by a Gaussian potential and 
denoted by Ufw(z'), are gradually inserted at the bound¬ 
aries of the simulation box at z = 0 and z = Lz with 
periodic boundary conditions maintained along all three 
axes. Since the structureless flat wall has a very short 
range ( 0 . 001 ( 7 ), few crystalline particles interact with it. 
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FIG. 2. (Color online) Crystal (blue) in contact with several 
layers of frozen-in crystalline particles (red), forming walls on 
either side. This snapshot represents an instataneons state 
of the system at equilibrium corresponding to the modified 
Lennard-Jones potential (mLJ) at p = 1.1 and T = 1.0. 


responds to the crystal in contact only with a flat wall, 
while at A = 1 the crystal is in contact with the frozen 
wall in presence of the short-ranged flat wall. 

The free-energy difference in each step is computed 
from the relation 

where the integration over A is computed using Simp¬ 
son’s rule. Finally, the interfacial free-energy is directly 
obtained from AFi and AF 2 via 7 = (Afi -b AF 2 )/A. 


As a result, the contribution of this step to the total in¬ 
terfacial free energy is negligible and does not affect the 
structure of the crystal near the walls. To integrate the 
equations of motion in molecular dynamics in presence 
of such extremely short ranged walls, a very small time- 
step is needed. This reduces the computational efficiency 
since the other forces in the simulation such as those be¬ 
tween the crystalline particles are more comparatively 
long-ranged (~ a ). However, by the use of a multiple 
time-step scheme , the computational overhead is 

reduced significantly and our simulations are slightly less 
than two times slower than those carried out with a single 
large time-step. 

In the second step, interactions of particles through 
the periodic boundaries along the 2 : direction are grad¬ 
ually switched off while interactions between the crystal 
and the frozen crystalline wall are gradually switched on, 
in presence of the flat walls. The final state is a crystal 
with periodic boundaries along the x and y axes and con¬ 
fined by frozen crystalline walls and the flat walls along 
the z direction. Since the free-energy difference per unit 
area corresponding to the first step is negligible the total 
contribution to 7 comes only from the second step. 

In the TI scheme adopted in this work, similar to our 
earlier works, the transformation from one state to the 
next is brought about by changing a switching parame¬ 
ter A, which couples to the interaction potential between 
the walls and the crystalline particles. The TI scheme 
consists of two steps. In step I, the A dependent Hamil¬ 
tonian has the form i?i(A) = where Uiv,{z) = 

aexp(—[(z — Zw)/1>]^), with a = 25.0fcBT and b — O.OOlcr. 
By switching A from 0 to 1, the system is transformed 
from a pure bulk crystal (A = 0, no wall), to a crystal 
in contact with an extremely short-ranged Gaussian flat 
wall (A = 1). In step 2, the A-dependent part of the 
Hamiltonian has different forms for the different inter¬ 
action potentials. For the LJ and inverse-twelfth power 
potentials, i? 2 (A) = (1 —A)^M*(r)-|-A^itu,(r), while for the 
hard-sphere potential, H2{X) = (1— A)®u*(r)-|-A^®®Mu,(r). 
Here, u*{r) represents interactions between the crys¬ 
talline particles through the periodic boundaries, while 
Uw denotes the interaction between a crystal and a frozen 
wall particle and is identical to the interaction between 
two crystalline particles. In the second step, A = 0 cor- 


V. SIMULATION 

To integrate the equations of motion, the velocity form 
of the Verlet algorithm [ 2 ^, [2^ was implemented with a 
multiple time-step scheme. For the mLJ and fsil2 poten¬ 
tials a smaller time-step Atsmaii = 0.00025r was used 
along with a larger time-step Atiarge = 0.004r, while 
for the hard-sphere potential, Atsmaii = 0.00025t and 
Atiarge = O.OOOSt were chosen. To maintain constant 
temperature, T = 1.0, every 200 time steps the velocity 
of the particles was drawn from a Maxwell-Boltzmann 
distribution at the desired temperature. 

We carried out Molecular Dynamics simulations in the 
NVT ensemble at various densities of the crystal at the 
temperature T = 1.0. Initially, an ideal fee crystal was 
put in the simulation box and it was ensured that the 
center of mass of the crystal along the z direction was 
at — Lz/ 2 , i.e. the two outermost crystalline layers are at 



FIG. 3. (Color online) Density profile of crystal in contact 
with a wall comprising frozen-in crystalline layers (solid line) 
and a wall consisting of particles arranged in an ideal fee lat¬ 
tice structure having the same density as the crystal (dotted 
line). The density profiles were averaged over 1000 indepen¬ 
dent configurations of the crystal in contact with the same 
frozen wall. The vertical dashed lines indicate the positions 
of the short-ranged flat walls at z = 0 and L^. The profiles 
were obtained for the mLJ potential at the density p = 1.1 
and the temperature T = 1.0. 
















































































































4 


the same distance from the boundaries at the two ends. 
At all densities, the crystal was first equilibrated at the 
desired temperature in order to chose the reference con¬ 
figurations for the pinned layers of crystalline particles 
forming the wall. 

The frozen wall is constructed by choosing particle po¬ 
sitions from crystalline slabs of width 2.5a (l.Ocr in case of 
the hard-sphere potential) from the boundaries at z = 0 
and z = Lz, at the two ends of the simulation box, 
from an equilibrium configuration. When equilibrating 
the crystal, every time the thermostat was imposed the 
average center of mass momentum of the crystal was set 
to zero. New particles fixed at these instantaneous equi¬ 
librated positions were then juxtaposed on the right and 
left sides of the simulation cell from the slabs at the left 
and right ends of width 2.5a or a [equal to the cut-off 
range of the interaction potential u{r)], by shifting their 
z-positions by —L^ and L^, respectively. 

The system sizes corresponding to the mLJ and fsil2 
interaction potentials consisted of iV = 6480 particles 
while for the hard-sphere potential, N = 3600, were cho¬ 
sen. Simulations carried out at larger system sizes yielded 
values of 7 in agreement with the smaller systems within 
the statistical errors. For each interaction potential, 7 
was computed at various crystal densities ranging from 
/9 = 1.0 — 2.0. At each density, 7 was calculated for sev¬ 
eral realizations of wall configurations. However, data 
corresponding to the different realizations were in close 
agreement with each other, within the numerical errors. 
In all values of 7 reported in the next section, the error 
bars are less than the symbol size and hence they are not 
specified. 


VI. RESULTS 

Figure [5] shows the simulation set-up of the system 
with the crystal in contact with several frozen-in layers 
of crystalline particles forming walls on both sides. It is 
clear from Fig.|2]that the structure of the crystal near the 
wall is identical to that in the bulk. This is also clear from 
the density profile shown in Fig. |3l There is no change 
in the density profiles in the layers in contact with the 
wall and those far away from it. We have also observed 
that the energy profile of the crystalline layers near the 
wall is the same as in the bulk indicating that the wall 
particles do not change the structure and energy of the 
crystal layers in contact with it. If, instead, the wall is 
constituted by particles attached to ideal fee lattice sites, 
the peaks of the density profiles near the wall are higher 
than those in layers away from the walls (cf. Fig. [3]), 
indicating that the structure of the crystalline layers in 
the vicinity of such an ideal wall is different from the 
bulk. 

The thermodynamic integrands corresponding to step 
2 of the TI scheme and pertaining to the different densi¬ 
ties are plotted as a function of A in Fig. 01 for the three 
interaction potentials. The contribution of the Gaussian 



FIG. 4. (Color online) Thermodynamic integrand as a func¬ 
tion of A for the second step of the TI scheme at various 
densities p for (a) the mLJ potential, (b) the fsil2 potential 
and (c) the 1256 potential. 



FIG. 5. (Golor online) Interfacial free energy vs. density for 
the (100) orientation of the crystal in contact with a frozen 
wall for the three interaction potentials. The dashed horizon¬ 
tal line indicates the zero axis. 
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flat walls in the first step is negligible (less than 0 . 1 % of 
7 ) and hence not reported. A negligible value of AFi 
indicates that the Gaussian flat walls effectively do not 
perturb the thermodynamics of the system. The contri¬ 
bution to 7 comes almost entirely from the second step. 
As Fig. |4] shows, the thermodynamic integrands corre¬ 
sponding to the two relatively long-ranged interaction po¬ 
tentials are similar, while that of the short-ranged hard- 
sphere potential is qualitatively different. The smooth¬ 
ness of the integrands leads to an accurate determination 
of the free-energy difference from numerical integration. 

In Fig. O we show 7 for the three interaction potentials 
as a function of the density. For the hard-sphere system, 
the maximum density we consider is p = 1.4 since the 
packing fraction at this density is ^ = 0.733, which is 
very close to the maximum theoretical packing fraction 
'('max ~ 0.74. For the other two soft sphere potentials, 
we are able to consider larger densities, up to p = 2 . 0 . 

Figure [ 5 ] shows a non-zero value for 7 as a function of 
p for all the three interaction potentials. The Helmholtz 
free-energy F is given hy F = U — TS, where U is the 
energy of the system, T the temperature and S the en¬ 
tropy. Since the potential energy of the system in contact 
with the wall is the same as in the bulk and the system 
is maintained at the same temperature throughout the 
thermodynamic transformation from the initial to the fi¬ 
nal states, a non-zero free-energy difference can only be 
attributed to a change in the entropy of the system. A 
non-zero free energy difference indicates that the frozen- 
in layer of crystalline particles imposes an external field 
on the crystalline particle, perturbing the thermodynam¬ 
ics as compared to the bulk. 

For the purely repulsive inverse-twelfth power potential 
and the hard-sphere potential, the interfacial free energy 
is always positive from the lowest to the highest densities 
considered. As can be seen in Fig. [U for rja > 1.0, the 
inverse twelfth-power potential is more repulsive as com¬ 
pared to the hard sphere interaction, resulting in a larger 
interfacial free energy. The positive excess free energy 
corresponding to both the purely repulsive potentials at 
all densities and the modified LJ potential at high den¬ 
sities originates due to a decrease of the entropy of the 
particles near the wall as they have less number of ways 
to arrange themselves near the wall and the thermal os¬ 
cillations of the particles around their mean positions is 
also suppressed in the vicinity of the wall. 

In case of the modified LJ potential, the interfacial 
free energy is positive only at high densities. As we go to 
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VII. CONCLUSION 
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